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The imaginary time path integral formalism is applied to a nonlinear Hamiltonian for a short 
fragment of heterogeneous DNA with a stabilizing solvent interaction term. Torsional effects are 
modeled by a twist angle between neighboring base pairs stacked along the molecule backbone. 
The base pair displacements are described by an ensemble of temperature dependent paths thus 
incorporating those fluctuational effects which shape the multisteps thermal denaturation. By sum- 
ming over ~ 10'^ - 10** base pair paths, a large number of double helix configurations is taken into 
account consistently with the physical requirements of the model potential. The partition function 
is computed as a function of the twist. It is found that the equilibrium twist angle, peculiar of 
B-DNA at room temperature, yields the stablest helicoidal geometry against thermal disruption of 
the base pair hydrogen bonds. This result is corroborated by the computation of thermodynamical 
properties such as fractions of open base pairs and specific heat. 

PACS numbers: 87.14.gk, 87.15.A-, 87.15. Zg, 05.10.-a 



I. Introduction 

Nonlinearities in DNA dynamics were first emphasized 
by Englander et al. [1] who interpreted the formation 
of temporary open segments of base pairs as mobile de- 
fects coherently propagating along the molecule back- 
bone. That seminal paper put forward the idea that 
the open configuration, allowing for hydrogen exchange 
with the solvent, involved torsional oscillations of the 
base pairs around the backbone axis. Hence the ther- 
mally activated excitations of the double helix could be 
twisted solitons spread over a length of about ten base 
pairs. Successive works [5] suggested that anharmonicity 
in hydrogen bond stretching modes could lead to self- 
trapping effects |3] and solitary-wave energy concentra- 
tion thus providing the energy flow required in RNA tran- 
scription [3]. Thus, soliton-like excitations were proposed 
to carry energy along the double helix in a way similar to 
the transport of biological energy occurring in a— helical 
protein molecules |5]- 

Later on, it was pointed out [6l [7] that DNA biological 
functioning does not necessarily requires energy trans- 
port whereas nonlinearities should be rather considered 
in the framework of models able to overcome the contin- 
uum approximation thus accounting for the discreteness 
of the structure PflUI- 

A remarkable amount of studies published over the last 
decades [Hi [12] has shed some light on the dynamics of 
DNA and the formation of fluctuational openings, the 
bubbles, which affect the thermal properties of the dou- 
ble helix including eventually its melting at high tem- 
perature. To tackle these issues, two classes of theoreti- 
cal methods have been developed. The first is based on 
the Poland-Scheraga model [131 [H] treating denaturing 
DNA essentially as a two state Ising-like chain of base 
pairs with regions of variable size, denaturated loops, 
opening temporarily due to thermal fluctuations. The 



second class assumes an Hamiltonian approach gener- 
ally based on the Peyrard-Bishop model 6^ in which the 
transverse stretchings between complementary base pairs 
are represented by a one dimensional, continuous vari- 
able bringing the advantage that also intermediate states 
in the DNA dynamics can be described. The Hamil- 
tonian approach, describing the system at the level of 
the base pairs, seems particularly convenient in hetero- 
geneous systems where the relative content of adenine- 
thymine{AT) versus guanine- cy to sine{GG) pairs |15j and 
sequence length (16| are key to determine thermodynam- 
ics and melting of the double helix. 

However, inside both theoretical approaches, different 
interpretations persist regarding the classification of the 
denaturation with some investigators finding a smooth 
crossover [TTj - EO] whereas others point to a sharp tran- 
sition driven either by self-avoidance effects [2TH23] or 
by nonlinear stacking along the molecule backbone [H] 
or by the finite range of the stacking itself [2B]. Denat- 
uration in synthetic homopolymer DNA displays a sin- 
gle melting temperature while heterogeneous DNA frag- 
ments denaturate in multisteps between about 310 and 
AlOK depending on the sequence and salt concentration 
of the solvent [261127]. 

Microscopically, experiments based on proton- 
deuterium exchange methods .28] have shown that 
the base pair lifetimes are of order of milliseconds at 
T ^ 305K with AT-pairs lifetimes being three times 
shorter than those of GC-pairs. Fluorescence correlation 
spectroscopy [29] has further revealed that the breathing 
modes have a time scale in the 20 — 100 microsec- 
onds range and bubbles of 2 to 10 base pairs open at 
T ^ 310K under low salt conditions. Base pair openings 
are thus highly localized at biological temperatures and 
remain such also in homopolymers. 

Theoretical modeling for DNA has then to incorpo- 
rate nonlinear dynamics and large fluctuational effects 
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which lead to base pair disruption and formation of bub- 
bles whose size varies with temperature and environment 
pH values [501 ISI]- While fully atomistic descriptions 
become computationally very heavy even for short frag- 
ments, mesoscopic models have proven capable to cap- 
ture the essentials of the interactions in DNA molecules. 
As fluctuations are expected to be strong in finite size 
sequences, I have recently proposed to apply the path in- 
tegral formalism [32 [33] to a one dimensional model, the 
Dauxois-Peyrard-Bishop (DPB) Hamiltonian [7], which 
incorporates nonlincarity also in the stacking potential. 
The idea underlying the method is that the base pair 
elongations with respect to the ground state are time 
dependent fluctuating paths which add to the partition 
function and cooperatively shape the denaturation tran- 
sition in flnite fragments. 

Including in the computation a great number of 
molecule configurations, I have found that the denatu- 
ration is a smooth crossover both in homopolymers and 
heterogeneous systems. The latter also display the well 
known multistep melting features in the specific heat 
plots mirroring the fact that segments of the molecule 
open at different temperatures with AT-rich regions driv- 
ing the bubble formation at relatively lower tempera- 
tures. Interestingly a recent study [34], based on the 
extended transfer matrix method [35] , which focusses on 
the same heterogeneous sequences considered in Ref . [53] , 
finds similar results regarding the temperature location 
of the main specific heat peaks. 

However the Hamiltonian model proposed so far suf- 
fers of a serious shortcoming as it does not account for 
the helicoidal structure which would have the immediate 
effect of bringing closer to each other non-consecutive 
bases along the molecule backbone [55J [37] . As a first 
step to describe helicity in the path integral method one 
may generalize the stacking Hamiltonian by introducing 
the angle of rotation between a base pair and the pre- 
vious one. In B-DNA at room temperature, one turn of 
the helix hosts about ten base pairs [35]. Accordingly 
the equilibrium twist angle, Oeq = 27r/10.4, is expected 
to provide the energetically most favorable configuration, 
the stablest one against thermal disruption of the base 
pair bonds. 

This Ansatz is a relevant benchmark for the compu- 
tational method based on path integrals and should be 
viewed as a constraint for DNA theories. In this work 
I develop the formalism for a generalized DPB Hamilto- 
nian in which the twist angle between neighboring bases 
appears as a free parameter. The Hamiltonian model 
also includes a solvent interaction term which realisti- 
cally simulates the formation of hydrogen bonds with 
the solvent thus stabilizing the denaturated state [39] . 
The partition function has been computed by varying 
the strength of the solvent factors. In Section II the 
model Hamiltonian is presented while the main features 
of the path integral formalism are outlined in Section HI. 
The thermodynamic results of the work are reported in 
Section IV where specifically the fractions of open base 
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FIG. 1: (Color online) (a) Schematic fixed planes picture for 
the right-handed helicoidal model. The blue filled circles de- 
note the pointlike base pairs stacked along the molecule axis 
with positive twist angle Q. The radial coordinate y„ de- 
scribes the pair mates displacement in the n — th base from 
the ground state. The dashed vertical axis corresponds to the 
y„ = configuration, the minimum for the one-coordinate 
potential VMiVn) + VsoiiVn) in Eq. |T|. There is no tilting 
between the planes, (b) The helix plane seen from above. 



pairs, the specific heat and the base pair paths ensem- 
ble are calculated as a function of the twist angle. Some 
conclusions are drawn in Section V. 



II. Hamiltonian Model 

The DPB Hamiltonian [7] for a system of N base pairs 
assumes the pair mates separation j/„ (for the n-th base 
pair) with respect to the ground state position as the 
relevant degree of freedom. The longitudinal base pairs 
displacements along the molecule backbone are neglected 
as they are much smaller than the transverse stretch- 
ings y„, hence the model Hamiltonian is essentially one- 
dimensional. The DPB Hamiltonian incorporates nonlin- 
earities both in the inter-base pair interactions modeled 
by a Morse potential and in the coupling between neigh- 
boring bases along the two strands. As these are viewed 
as two parallel chains, the DPB model does not account 
for helicity. 

Here, see Fig. [T] I introduce a twist angle 9 between 
adjacent bases, n and n — 1, along the DNA backbone 
thus generalizing the DPB Hamiltonian to the following 
expression 
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N r -9 



n=l 



+ Vs{yn,yn~l) + VMiVn) + VsoliVn) 



K 



Vs{Vn,yn-i) = y5«,n-i(y« - 2y„y„- 1 COS 6* + 1 ) 

5«.n-l = 1 + pexp[-Q;(j/„ + Vn-l)] 

VMiVn) = Aj(exp(-a„j/„) - l)^ 
VsoiiVn) = -D,Js{ta.nh{yn/ls) - l) • 

(1) 

H is the base pair reduced mass which is related to 
the harmonic stacking coupUng K hy K = fii/^, v be- 
ing the frequency of the phonon mode. As the DPB 
model is homogeneous with regard to the stacking inter- 
actions described by the potential Vs(?/„, /i takes 
the same value both for GC- and AT- bases. Also the 
anharmonic (positive) parameters p and a are assumed 
independent of the type of base at the n and n — 1 sites 
and I have checked that such homogeneity assumptions 
don't change substantially the thermodynamical proper- 
ties in the considered denaturation range. Instead, the 
latter are essentially determined by the size of K [ID] 
as it will be shown in Section IV. A. Hereafter the values 

o — 1 n — 2 

p— 1, a = 0.35A and if = QOmeVA are taken con- 
sistently with previous investigations on the DPB Hamil- 
tonian [331 HI] . Different values for the stacking pa- 
rameters are also found in the literature [42] . 

Note that p and a have a special role in driving the co- 
operative behavior of the system towards denaturation: 
when the molecule is closed, y„ , yn-^i ^ for all n, 
and the effective coupling is K{1 + p). If, however, a 
fluctuation causes either ?/„ > a^^ or > a^^, then 
the hydrogen bond between pair mates loosens and the 
base moves out of the strand axis. Accordingly the tt elec- 
trons overlap in the base plateaus is reduced, the binding 
between neighboring bases along the strand weakens and 
the effective coupling drops to K (in Eq. ([I]), .g„,„_i ~ 1). 
As a consequence, also the next base tends to move out 
of the stack thus propagating the fluctuational opening. 
This explains the link between anharmonicity and coop- 
erativity leading to bubble formation and eventually to 
denaturation in the DPB model. 

The relevant feature is the torsion, due to the cos 61 
term, in the stacking potential Vs{ymyn~i) in Eq- 
While such torsional effects have been included in sev- 
eral mesoscopic studies, ranging from molecular dynam- 
ics simulations [35] to extensions of DPB model [U] 
and of Ising-like model [15] , the role of the helicoidal ge- 
ometry on the denaturation patterns is still unclear. To 
shed light on this issue I use the path integral method 
to study the thermodynamics of the system in Eq. ([T]) 
as a function of 9. Specifically a set of double helices is 
considered, each with TV base pairs whose arrangement 
along the strands is determined by 9. The latter is taken 
as an input parameter. Only positive 9 are considered 
to avoid non helicoidal structures with alternating ±9 



between consecutive base pairs along the molecule back- 
bone: these zig-zag structures may also minimize the en- 
ergy as Vs{yn,yn-i) is even in 9. 

Strictly speaking, the stacking varies with the torsion 
of the molecule thus the couplings should also bear a de- 
pendence on 9, an effect which has not been quantified yet 
and it is here neglected. It is emphasized that the model 
m Eq. is stiU one-dimensional as the rotational de- 
gree of freedom is site independent. Accordingly also the 
fluctuations around the equilibrium angle in the coiled 
structure are not accounted for. Moreover, as there is no 
tilting for the base pair planes in Fig. [T] bending effects 
P5] are not included in the model. 

The base pair hydrogen bonds are modeled by the 
Morse potential in Eq. ([I]): Va/(?/„) incorporates hetero- 
geneity in the sequence through the site dependent pair 
dissociation energy Z)„ and the inverse length a„ which 
sets the potential range. As the energy per hydrogen 
bond is about 15meV, I set for AT— and GC— bases 
respectively Dat — 30meV and Dqc — AbmeV. AT- 
base pairs have larger displacements than GC-base pairs, 
then the inverse lengths are taken as qat = 4.2A and 
o-GC — 5^ while slightly different values can be found 
in the literature [47j|48]. Note that the plateau in Va/(2/„) 
has a relevant physical implication: once all base pair dis- 
placements are larger than a~^, the open strands can go 
infinitely apart with no energy cost. This means that 
the Morse potential does not account for strand recom- 
bination events which instead occur in solutions [5U1 25] 
and whose rate depends on the proton concentration in 
the solvent. This drawback has been circumvented by 
several techniques which aim to restrict the configura- 
tion space thus keeping the transverse stretchings within 
a finite range [3S]. While the path integral method nat- 
urally operates a truncation of the path configuration 
space [33] , I treat here the problem analytically by adding 
in Eq. ( 1 ) a solvent interaction term Vsoi (Un) (as proposed 
in Refs. |39Ll5n] ) modulated by a barrier factor fs and by 
a length Is setting the range of the potential. 

The solvent term has the effect to enhance (by fsDn) 
both the energy of the equilibrium configuration and the 
height of the barrier below which the base pair is closed. 
It is known that the melting temperatures depend log- 
arithmically on the sodium concentration [iVa"''] in the 
solvent [51] [55] . As the melting temperatures scale essen- 
tially linearly with the Morse potential barrier, also the 
latter can be assumed to vary logarithmically with [A^a^] 
[16) . This allows us to establish empirically a quantita- 
tive link between fs and [iVa+]. For instance, taking 
fs = 0.1, one gets [Na+] ^ 0AM. Instead, the length Is 
is taken so as to sample the effect of the solvent potential 
on the ensemble of the path displacements determined in 
the next Section. 

In Fig. [2] Vm (j/„) is plotted together with the super- 
imposed effect due to Vsoiiyn)- The potential parameters 
are those for GC-base pairs and also the bare Morse po- 
tential is shown. As a main feature VMiyn)+ Vsoiiyn) dis- 
plays a hump whose maximum determines the threshold 
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FIG. 2; (Color online) Sum of Morse Vm and solvent Vsoi 
potentials versus guanine- cytosine base pair separation for 
two values, in (a) and (b) respectively, of the solvent barrier 
factor fs. Is (in A) tunes the width of the solvent barrier. 
Also the bare Morse potential is plotted: the parameters D„ 
and Un are taken for GC-base pairs. 



stretching around which a base pair may first temporarily 
open and then either re-close or fully dissociate. Looking 
at the case Ig = 3 A in Fig. [2|b), the maximum occurs 
at ~ 0.8 A. Accordingly, in the range 0.8 < j/n < Is, the 
pair mates are pulled away from each other and, for base 
pair stretchings larger than 1^, the hydrogen bond with 
the solvent is established. In this sense the solvent inter- 
action stabilizes the inter-strand open configurations. 



the statistical partition function, the integration being a 
trace integration. In the calculations, not all histories 
can be accounted for and, given the specific problem, 
one has to select the suitable class of paths which mainly 
contribute to the physical properties. 

I have applied the path integral method to the discrete 
model in Eq. ([T]) by introducing the idea that the N base 
pair displacements yn can be described by one dimen- 
sional paths x{Ti), the latter being periodic functions of 
the imaginary time x^, x{Ti) — xiji -\- (3). The index i 
numbers the base pairs along the r-axis. In fact, there 
are -I- 1 base pairs in Eq. ([T]) but the presence of an 
extra base pair is remedied by taking periodic bound- 
ary conditions, yo = UN, which close the finite chain 
into a loop. This condition is incorporated in the path 
integral description as the path is a closed trajectory, 
x(0) — x{j3). Hence a molecule configuration is given by 
A^^ = A" paths and, in the discrete imaginary time lat- 
tice, the separation between nearest neighbors base pairs 
is At = P/Nr- Then, Eq. ([I]) is mapped on the time axis 
as follows: 



Vn x{Ti) 

yn-i x{t') 
t' At 
n= l,..,N 
i= l,..,Nr- 
Ti = ; 



13. 



(2) 



Accordingly also the real time derivative y„ maps onto 
the imaginary time derivative x{t) as: 



III. Computational Method 

In the path integral formalism, the zero temperature 
evolution amplitude between two points, say "a" and "b" , 
is a sum over all histories along which a system can evolve 
in going from "a" to "b" during the time t. Each history 
is weighed by a phase factor, the exponential of the ac- 
tion associated to a given path [S3]. The formalism is 
extended to the finite temperature case, by performing 
an analytic continuation which defines the axis of the 
imaginary time r with t — it. As t € [0, /3] and /3 is the 
inverse temperature, the imaginary time path integral 
permits to derive the thermal properties of the system 
by weighing the contributions to the action due to the 
particle paths x{t) [51] . 

Accordingly the statistical partition function, which 
can be viewed as an analytic continuation of the quan- 
tum mechanical partition function, is given by an integral 
in the path phase space. Each path is weighed by a prob- 
ability factor exp{—A[x{T)]) where the Euclidean action 
A[a;(r)] replaces the mechanical canonical action of the 
zero temperature case. Only closed paths contribute to 



dyn 
dt 



dx 
d^ 



Eq. ([3]) is consistent with the replacement 



(3) 



(4) 



which is justified in the classical regime appropriate 
to DNA denaturation. Eq. (|4| is also used to solve the 
pseudo-Schrodingcr equation for a Morse potential [55] 
obtained from Eq. ([I]) in the large K limit (with p = 
and fs =0)m- 

Due to their periodicity property the paths can be ex- 
panded in Fourier series with cutoff Mp 



Mf 

x{Ti) = xo + ^ \a,n cos{uJmTi) + h,n sin(w„Ti) 

m— 1 

— 2m7r/^ 



(5) 



and this introduces the following physical picture: 
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a) given a set of coefficients {xq, am, &m}, tlie N base 
pairs are represented by the configuration {x{Ti), i = 
1,.., N^). 

h) A set of coefficients corresponds to a point in 
the path configuration space thus, samphng the latter 
amounts to build an ensemble of distinct configurations 
for the system. As this is done for any temperature, we 
have a tool to describe the base pair thermal fluctuations 
around the equilibrium {x{Ti) ^ 0). 

c) In principle the configurations ensemble for the 
DNA fragment is infinite as it may include any possible 
combination of Fourier coefhcients. For practical pur- 
poses some physical criteria intervene to select computa- 
tionally the path coefficients defining a molecule config- 
uration and contributing to the partition function. This 
poses a restriction on the ensemble size. 

Such criteria are of two types: First, the hydrogen 
bond potential in Fig. [2] excludes too negative base pair 
stretchings due to the hard core which mimics the repul- 
sion between negatively charged sugar-phosphate groups. 
Also too large displacements, x{Ti) ^ Ig, are eliminated 
as they have no effect on the free energy derivatives. 
Once the pair mates are unbound and tied to the sol- 
vent molecule there is no physical reason to further pull 
them to any particular direction. Second, the ensemble 
of paths has to be consistent with the thermodynamics 
laws. This means that the numerical code selects, at any 
temperature, a path ensemble and evaluates the entropy 
of the DNA fragment. If the entropy is growing versus 
T, the code proceeds to the next temperature step oth- 
erwise a new partition is performed in the Fourier coef- 
ficients integration, a new path ensemble is selected and 
the entropy is recalculated. This is done at any T until 
the macroscopic constraint of the second law of thermo- 
dynamics is fulfilled throughout the whole investigated 
temperature range. I emphasize that the method does 
not put any constraint on the shape of the entropy ver- 
sus T- plot aside from the requirement that the entropy 
derivative has to be positive. 

It follows that the path ensemble is a dynamical object 
accounting for the manyfold of molecule configurations 
which enter the thermodynamical calculation. The size 
of the ensemble is a measure of the cooperativity degree 
of the system. By increasing T, some base pairs may 
open and cooperatively lead to bubble formation along 
segments of the double helix. Accordingly the ensemble 
size is expected to grow versus T. I define Nf,f f the num- 
ber of sets of Fourier coefficients in Eq. ([s]) , selected at a 
given temperature, which fulfill the criteria above. Then 
Neff is the number of possible trajectories for the i — th 
base pair while the ensemble size for the whole fragment 
is measured by Nt x N^ff. This is a key parameter in 
our analysis as it will be shown below. 

So far, one molecule which may exist in N^f f configura- 
tions has been assumed. In a picture closer to the experi- 
mental viewpoint, one may consider an ensemble of iden- 
tical molecules (with N^. base pairs) . Each molecule may 
exist, at a given T, in a configuration specified by a set of 



base pair displacements hence by one point {xq, am, bm} 
in the path configuration space. This establishes a biu- 
nivocal correspondence between the molecules ensemble 
and Nef / while the overall base pair ensemble size is mea- 
sured by N-r X Nef / . 

Applying the mapping technique in Eqs. ([2]), (|3| to 
Eq. ([T]), the classical partition function for the DNA 
molecule in the solvent is written as 

Zc = (j) T)xexp[-l3Ac{x}] 
Ac{x)= 5][|i(r,f + T/5Wr,),a:(r')) + 



VM{x{Ti))+Vsol{x{n)) 




(6) 

where is the thermal wavelength which, by virtue of 
Eq. Q, takes the expression A^ = ^Ji: / (3K. 2)x is the 
measure of integration and § indicates that the paths 
x[t) are closed trajectories [551 [57]. 

IV. Thermodynamical Results 

The theory is applied to a fragment with Nr ~ 100 
base pairs whose sequence is: 

GC + QAT + GC + HAT + SGC + AT + AGC + 
AT + AGC + AT + 8GC + [49 - 100] . 

(7) 

Due to the prevalence of AT- base pairs, bubbles may 
open in the leftmost and, more likely, in the right part. 
The left portion of 48 base pairs has the same sequence 
as the L48AS fragment of Ref. [HHj although our model 
cannot distinguish, for instance, a AT- followed by AT- 
along the backbone from AT- followed by TA-pair. Such 
differences may shift considerably the melting tempera- 
ture in short fragments [SHI HO] with effects which are 
quantitatively not fully understood. Coherently with the 
notation used in Ref.[33 , the sequence in Eq. ([7| will be 
hereafter named L48AT22. 



A. Fractions of Open Base Pairs 

The melting temperature is experimentally defined as 
the temperature at which half of the molecules in the 
sample are in the double-helical state and half are in the 
single strand, random-coil state. As explained in the pre- 
vious Section, this amounts to say that half of the con- 
figurations Neff for the fragment in Eq. ^ are closed 
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and half are open. When base pairs dissociate the UV 
signal changes quite abruptly ^Glj. However, once the 
UV signal measures that half of the base pairs are open, 
this may indicate either that half of the molecules in the 
sample are open and half are closed or that all molecules 
arc half-open. Accordingly absorption methods cannot 
distinguish intermediate states for a single molecule con- 
figuration which, instead, would be quite interesting in 
order to understand the nature of the melting transition. 
New techniques based on quenching of single strands are 
becoming available to trap intermediate states [55] . 

On the theoretical side, the problem arises to define 
when a configuration is open or closed. There is necessar- 
ily some arbitrariness intrinsic to the Hamiltonian model 
as this is expressed in terms of base pair stretchings which 
are continuous variables. Accordingly one may assume 
that a configuration is open when all Nr base pairs are 
larger than a certain threshold which however cannot 
be set univocally. For these reasons even large discrepan- 
cies regarding the choice of C are found in the literature 
[621163]. 

To tackle these questions I have computed the fraction 
of open base pairs F^p for the system in Eq. ([tI taking 
different C and searching for a range of values which yield 
a description of the denaturation qualitatively consistent 
with that provided by thermodynamical indicators such 
as the specific heat fBT . This reasoning is inspired by the 
observation that the fraction of closed base pairs, 1 — Fop, 
is a measure of the system internal energy hence dFop/dT 
is proportional to the specific heat. This holds for a ho- 
mogeneous chain but also in heterogeneous DNA the spe- 
cific heat is an indicator of the melting as it displays sharp 
peaks at the temperatures where various parts of the se- 
quence open. Thus some correlation is expected between 
Fop and specific heat plots versus temperature. 

Using Eq. Q the fraction of open base pairs is given, 
in terms of path integrals, by 

1 

Pop - ]^E^(<^(^^)>-0 

^ i—l 

< x{t^) >^ Z^^ (j) Dxx{T^)exp[-(3Ac{x}] . (8) 

The Heaviside function is justified by the sharp 
increases in the UV signal at the denaturation steps. In 
fact, in Eq. ([8|, x{Ti) is averaged over the configuration 
ensemble and only < x{Ti) > is directly confronted to 

Then, to be rigorous, every configuration might have 
at least a base pair such that, x(ri) < Ci (and therefore 
to be classified as closed) whereas the average configura- 
tion might still be open as, < x{Ti) > — ^ > 0, for any 
i. With this caveat and aware that there exist alterna- 
tive definition for Fop [IHl US] directly relating the base 
pair stretching to the opening threshold, Eq. ^ is here- 
after used as it better captures the multiple steps in the 
denaturation. 

I have first evaluated the role of the solvent potential 




Temperature (K) Temperature (K) 



FIG. 3: (Color online) Fractions of open base pairs, calculated 
via Eq. Q, for two solvent barrier factors: (a) fs = 0.1, (b) 
fa — 0.3. la is in units A. Three thresholds ^ are assumed: 
(a) C = 0.8, 1.0, 1.5A, (b) C = 0.6, 0.8, l.OA. The plots 
show the fractions of base pairs larger (>) than versus tem- 
perature. 



in Eq. ([8]) assuming an untwisted DPB model, thereby 
9 = in Eq. ([T]) . The results are reported in Fig. [s] for 
those two barrier factors fs considered in Fig. [2| First 
one notices that fs has to be sufficiently large in order to 
make also the effect of Is appreciable. Fop is calculated 
taking three ( in Fig. |3]^a) and Fig. [sjjb) respectively. 
Two values, ^ = 0-8, 1^, are common to both figures. 
It is seen that the solvent has a stabilizing effect as, for 
the same ^ but larger fs, Fop attains the unity at higher 
temperatures. In Fig. [sjja), all averaged paths become 
larger than ( = 0.8A at T ~ 3Q0K whereas this event 
occurs, m Fig.[3];b), at T - 350K: a sizeable shift which 
points to the importance of the solvent in modeling DNA 
fragments with the purpose to fit experiments. Certainly 
the fact that all average base pair paths are larger than, 
say ( — 0.8A, at a given T does not necessarily mean that 
the sequence is melting at that T. However, had we to 
know the melting temperature for a specific fragment, 
we would able to estimate via Eq. ^ the C value such 
that Fop = 1 at that Tm ■ In this view the computational 
method can select a reliable threshold for the base pair 
elongations with respect to the equilibrium. Above such 
threshold the average configuration may be considered 
as open. In substantial agreement with Refs.fW, !M] our 
numerical work suggests that a range C ^ [0-5 — 1.0]A 
is appropriate to the purpose of modeling a multistep 
denaturation which occurs at Tm ^ [300 — 390]if . This 
range is chosen for the next calculations together with 
fs = 0.3 and Is = 3.0l. 

Next we turn to the core of the work by estimating 
the effect of the twist on Eq. ([s]). Beside the untwisted 
case {6 — 0), two angles 9 = 0.60707, 0.30353 rad are 
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FIG. 4: (Color online) Fractions of open base pairs versus 
temperature calculated according to Eq. ([8|. The green filled 
circles and stars indicate the fractions of base pairs larger than 
the thresholds ^ — 0.6 A and = 1.0 A respectively. Three 
twist angles (in rad) are assumed in the stacking potential 
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The backbone coupling is A" = GOmeVA 



Is 



considered which correspond to accommodate in one turn 
of the helix about 10 and 20 base pairs respectively. The 
results are plotted in Fig. [4] vifhere Fop for the sequence 
L48AT22 is calculated versus temperature assuming two 
threshold values. 

Fop varies strongly with 9 and gets larger in the un- 
twisted model suggesting that the latter may denatu- 
rate at lower temperatures. Setting ( = 0.6A one ob- 
serves that, by increasing 9, Fop = 1 is attained at 
T - 275, 300, ibm respectively This docs not allow 
us to conclude that the case 9eq = 0.60707 rad corre- 
sponds to the stablest configuration for the fragment but 
Fig. [4] suggests that twist angles smaller than 9f,q would 
procluce an helicoidal geometry which undergoes large 
fluctuational effects already below room temperature. 

So far the melting profiles have been computed taking 
a constant harmonic stacking K. As the overall stability 
conditions for the double helix [ID] are determined by the 
base stacking it is worth analyzing Fop as a function of K 
in the presence of the solvent potential [66] ■ The results 
are shown in Fig. [5] for the twists considered in Fig. [4| 
The threshold is set at C = O.SA which corresponds to the 
maximum in the barrier height for GC-pairs due to the 
solvent (see Fig. 2]). The strength of K has a large weight 
in the untwisted structure (Fig. [sja)): for instance, fo- 
cusing on T ^ 34Gi4r, a K enhancement from 60 to 

o —2 

^QmeVA produces a ~ 40 % drop in Fop, whereas even 
70 % of the average base pair displacements become 
smaller than C = 0.8A in the case K ^ IQQmeVA 
The effect of the backbone stiffness is much less dramatic 
in the twisted configurations as the twist on its own gives 
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FIG. 5: (Color online) Fractions of open base pairs, larger 
than the thresholds C, = O.SA, calculated through Eq. ((sk for 
different backbone couplings K (in units of meVA ). The 
three twist angles 6 of Fig. J4] are assumed in (a), (b) and (c) 
respectively, h is in units A. 
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stability to the system. This can be observed in Fig. 
and, more remarkably, in Fig. [sj^c). In the latter, at 
T ^ 340K, Fop drops by only ^ 10 % in going from 60 

n —2 

to IQOmeVA . Note however that the melting feature 
at T ~ 380K displayed by the K = 60 plot is smeared 
at larger K confirming that accurate determinations of 
the stacking parameters in heterogeneous sequences are 
fundamental for quantitative predictions of the melting 
profiles. The backbone coupling also affects the degree 
of cooperativity in the denaturation transition as N^f j is 
strongly reduced by increasing the stiffness via K [H5] . 

B. Specific Heat versus Twist 

Computation of Eq. ^ yields the free energy A = 
—/3~^ In Zc hence the thermodynamic properties of the 
fragment in Eq. ([7| . For any choice of the model param- 
eters the entropy is always found to be a continuously 
growing function of temperature confirming our previous 
analysis regarding the character of the denaturation: this 
is a smooth crossover which, in heterogeneous fragments, 
appears in niultisteps driven by the AT-rich regions of 
the sequence. While the entropy plots (similar to those 
in Ref. [33]) are not displayed here, I focus on the spe- 
cific heat which is shown in Fig. |6] both for the untwisted 
case and two twist angles, 9 ~ 0.60707, 1.57079 rad. The 
harmonic backbone coupling is set at if = GOmeVA 

The specific heat for 9 = 0.30353 rad (the value taken 
in Fig. |4]) is not reported as it would almost overlap with 
the untwisted curve. The plot with 9—0 shows two 
shoulder peaks at T ~ 275, 300K which are consistent 
with the sharp increase in Fop, for C — 0.6A, below room 
temperature. The main peak at T ^ 330A' signals the 
denaturation of a large portion in the ensemble configu- 
ration and, confronting with Fig.[3](b), the corresponding 
opening threshold is C = O.SA. However not all molecule 
configurations have opened at T ^ 330K: some more do 
it at T ^ 388K as shown by the further shoulder peak 
witnessing that the solvent has effectively stabilized the 
system up to this value. In fact the T ~ 388K- peak is 
smeared by switching off fs- Note that the two dips at 
T ^ 290, 380K don't have any thermodynamical mean- 
ing: both can be smeared by slightly enhancing the path 
ensemble size rather indicating an overall lower stability 
for the ladder configuration of the untwisted DPB model. 

Similar is the trend for the 9 = 1.57079 rad specific 
heat presenting small size peaks below room tempera- 
ture, a main peak at T ~ 336K and a shoulder peak at 
T ^ 362K. Instead, the 9eq- plot does not show any low 
T peak indicating that the system is substantially sta- 
ble up to room temperature. The main peak occurs at 
T ^ 320A' where ~ 65% of the base pair configurations 
have average values larger than C = 0.6 A. All configura- 
tions exceed such threshold at T ~ 350K but this does 
not suffice to achieve complete melting as some base pairs 
open at higher T once their stretchings exceed a larger 
( (~ 0.8 A). These findings are qualitatively consistent 
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FIG. 6; (Color online) Specific Heat versus temperature cal- 
culated via Eq. ([6| by varying the twist angle 9 (in rad). The 
backbone coupling is K — QOmeVA 

with a neutron scattering investigation of B-DNA f2Q] 
showing the persistence of closed base pairs regions well 
inside the denaturation regime while the melting transi- 
tion appears overall smooth. 

The specific heat curves reveal multistep melting with 
fluctuational openings for all investigated twist angles. 
However only the system with does not present de- 
naturation peaks below room temperature pointing to a 
higher stability for such helicity. This result, in itself 
peculiar, is a motivation to analyze further the model 
dependence on the twist. 



C. Base Pairs Path Ensemble versus Twist 

To complete the study, I consider the size of the con- 
figurations ensemble for the base pairs in the fragment 
which is defined by x Nf,ff. It says how many dis- 
tinct path configurations participate to Zc at any tem- 
perature thus providing a measure of the degree of co- 
operativity in the system. In homogeneous DNA, close 
to the melting transition, such degree grows sharply [32] 
as many base pairs are involved in the event at about 
the same T. In heterogeneous DNA the cooperativity 
degree grows continuously as the transition is spread out 
in niultisteps [57]. Although the T locations of the de- 
naturation steps is not precisely signalled by Nr x N^ff, 
the more stable is the heterogeneous fragment against 
denaturation the smaller N^- x Nefj is expected to be. In 
a way, this fundamental parameter for the path integral 
method measures the molecule resilience to the thermally 
induced perturbation. In Fig. [Tj A'^^ x A^e// is plotted 
versus T for the untwisted case and four twist angles. 
In all cases the calculation starts at T = 260Ar including 
~ 12 X 10® paths which ensure numerical convergence for 
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FIG. 7: (Color online) Total number of paths used in the 
computation of the DNA fragment thermodynamics for the 
untwisted case and four values of the twist angle (in rad). 
The backbone coupling is K = GOmeVA 



Zc- Then, Nr x N^fj is re- normalized upwards versus 
temperature according to the method described in Sec- 
tion III. As a main result, the helicity 9eq is indeed the 
value providing more stability to the fragment. For all 
other considered twists 6, Nr x Neff grows faster versus 
T indicating that the double strand molecule is more eas- 
ily disordered by thermal effects on the hydrogen bonds. 
This means that denaturation steps may occur at lower 
T than for the 0eq case although no proportionality re- 
lation between Nr x N^f / and twist angle can be gener- 
ally inferred. Note in fact that the case 9 = 1.57079 rad, 
which displays partial openings already at low T (Fig. |6]) , 
also shows the most rapid growth in the size of the path 
ensemble. To capture these features one has to perform 
summations over ~ 10® base pair paths in the upper por- 
tion of the considered T range, a highly time consuming 
computational task. 

The twist 9eq is found to yield the highest stability 
also for other sequences obtained by varying order and 
relative content of AT- and GC- pairs in Eq. ([7| . 

This study has sampled different twist configurations 
searching for the most energetically advantageous. The 
mesoscopic Hamiltonian in Eq. ([T]) does not account for 
torsionally constrained or enzyme driven supercoiling ef- 
fects [55]. However, by virtue of the closed boundary 
conditions on the path displacements, it is worth to in- 
vestigate whether the model may provide informations 
for the occurrence of supercoiling in small circular DNA 
molecules [69] . Say a the superhelical density which 
measures the degree of supercoiling with respect to the 
relaxed molecule: as natural DNA is negatively super- 
coiled, underwound configurations should be energeti- 
cally more favored than the overwound ones |70| . As- 
suming zero number of writhes, Nr x N^ff is computed 



by varying the number of twists Tw with respect to the 
equilibrium, {Tw)o = 10. This allows to monitor the co- 
operativity degree tuning the twist angle 9 around 9eq- 
For |Tw — {Tw)q\ — Z with integer Z, one gets two 
angles 0i=i,2- E.g., Z = 1 implies: a) Tw = 11 and 
6*1 = 0.69115 hence, a > 0; b) Tw ^ 9 and 92 = 0.56548 
hence, cr < 0. As A'^^ = 100 for the sequence L48AT22, 
the helical configurations obtained by taking Z € [1,3] 
will suffice to our purposes. 

With the potential parameters of Fig. [T) the model 
indeed accounts for the asymmetry between positive and 
negative supercoiling proving the energetic convenience 
of the (T < configurations. 

In fact, it is found that: 



(9) 



Z = 1 

N,ff{9,)^ N,ff{92)^ Neff{9eq) 
throughout the whole temperature range. Instead: 
Z = 2,3 

Neffi9l > 9,q) > N,ff{92 < 9,q) ^ N,ff{9,q) 

(10) 

in the upper portion of the temperature range. 

This amounts to say that small twist changes, i.e. 
Tw = 9, 11, do not perturb the thermodynamics of the 
relaxed molecule but once 9i= 1^2 differ significantly from 
6eq, i.e. for Tw = 8, 12, it is the underwound con- 
figuration to be naturally preferred: a moderate helix 
unwinding can sustain local denaturation bubbles which 
absorb the released twisting strain thus allowing for an 
overall molecule stability. Certainly, a too large unwind- 
ing (cr — —1) favors the complete strand separation and 
leads to melting. 

These findings are consistent with careful single 
molecule micromanipulation experiments |68| and the- 
oretical predictions based on the Poland-Scheraga model 
|71) . mostly in regard to the asymmetric system depen- 
dence on the sign of a. At this stage however our method, 
applied to the sequence in Eq. ([7|, does not display any 
relevant physical effect by varying the superhelical den- 
sity in the biologically interesting range, \a\ < 0.1. 



V. Concluding Remarks 

The thermodynamics of a short fragment of heteroge- 
neous DNA has been studied by the path integral method 
focusing on the temperature range in which denaturation 
takes place. The model Hamiltonian contains a solvent 
interaction term that enhances the base pair dissocia- 
tion energy and stabilizes the hydrogen bonds between 
complementary strands. Overcoming the limitations in- 
trinsic to previous investigations which model the dou- 
ble helix through two parallel chains, I have introduced 
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as a free parameter a twist angle between two adjacent 
base pair stacked along the molecule backbone. In B- 
DNA the strands twist, around the molecule axis, about 
once every ten base pairs hence, the equilibrium twist an- 
gle is 9eq — 0.60707 rad. Assuming an adenine-thymine 
rich sequence of one hundred base pairs, I have studied 
whether the path integral model can capture an energetic 
advantage for such equilibrium geometry. The base pair 
displacements are considered as imaginary time depen- 
dent paths and the ensemble of good paths which con- 
tribute to the system partition function has been built 
consistently with some physical constraints, markedly the 
shape of the hydrogen bond potential and the second law 
of thermodynamics. Such path ensemble is a tempera- 
ture dependent representation of the possible molecule 
configurations and it accounts for those fluctuational ef- 
fects which are key to the DNA dynamics, mainly in short 
fragments. In general, close to the denaturation, the path 
ensemble size grows as the system becomes more coopera- 
tive and more paths participate to the transition. Thus, 
the ensemble size measures the overall system stability 
against thermal disruptions of the hydrogen bonds. By 
varying the twist angle, it is found that the stability is 
higher precisely for the torsion defined by O^q and this 
conclusion does not depend on the specific sequence. The 
Ansatz made in the Introduction is then proved in the 
light of the path integral approach here developed. 

Two major physical effects appear in our model and 
compete on the energy scale: the base pair displacements 
should become large enough to yield those fluctuational 
openings which are peculiar of the double helix dynamics 
and, at the same time, too large base pair fluctuations 
should be discouraged as they destabilize the system al- 
ready at room temperature. The stacking potential has 
the capability to select the appropriate fluctuation am- 
plitudes as a function of the twist. 

In untwisted or small twist models, 9 <^ 9eq, even large 



fluctuations are possible as the corresponding stacking 
energies, Vs{yni ^n-i), are low on the energy scale set by 
the one coordinate potential VMiVn) + ^oz(2/n)- Hence, 
such large fluctuations do contribute to partition function 
and thermal properties. It follows that molecules with 
small twist have scarce stability and begin to denaturate 
below room temperature. 

On the other hand, in large twist models with 9 ^ 
9eq, even small fluctuations have stacking energies which 
are higher than the hydrogen bonds dissociation energy 
hence, their contributions to the partition function be- 
come vanishingly small. Accordingly large twist mod- 
els have scarce flexibility and do not account for those 
fluctuational openings that are vital to the molecule. In 
this view 9eq emerges from the path integral computa- 
tion as the energetically most convenient twist for the 
double strand configuration since it provides a right bal- 
ance among the competing tendencies of this complex 
system. 

For \9 — 9eq \ ^ 0.15, the unwound double helix displays 
an energetic advantage over the overwound one indicat- 
ing that the opening of local denaturation bubbles is a 
natural strategy to stabilize the system absorbing tor- 
sional strain. 

I have also computed the fractions of open base pairs 
with reasonable opening thresholds and the specific heat 
for various twist angles: all the obtained results provide 
as much independent as consistent indications that the 
helicoidal geometry of B-DNA has indeed a higher sta- 
bility, albeit maintaining that flexibility associated with 
the nonlinear character of the intra- and inter-strand in- 
teractions. Further investigations on these issues are due 
to come. In particular the twist should be treated as a 
rotational degree of freedom in a two dimensional path 
integral description thus incorporating also torsional fluc- 
tuation effects around the equilibrium geometry. 
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